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We present a polynomial Hybrid Monte Carlo (PHMC) algorithm as an exact simulation algorithm with dy- 
namical Kogut-Susskind fermions. The algorithm uses a Hermitian polynomial approximation for the fractional 
power of the KS fermion matrix. The systematic error from the polynomial approximation is removed by the 
Kennedy-Kuti noisy Metropolis test so that the algorithm becomes exact at a finite molecular dynamics step size. 
We performed numerical tests with Nf=2 case on several lattice sizes. We found that the PHMC algorithm works 
on a moderately large lattice of 16^ at (3=5.7, m=0.02 (mps/mv~0.69) with a reasonable computational time. 



1. Introduction 

The low energy QCD dynamics in the real 
world will be understood by lattice QCD with 
three-flavors of dynamical quarks. Several efforts 
have been spent to develop exact numerical algo- 
rithms with an odd-numbers of the Wilson type 
quarks [Q. 

The Kogut-Susskind (KS) fermion is an at- 
tractive formalism since the numerical simula- 
tion with much lighter quark masses are possi- 
ble thanks to the remnant chiral symmetry. Al- 
though lattice QCD with the two- or single-flavor 
KS fermions can be defined by taking the frac- 
tional power of the KS fermion, efficient exact 
algorithms are not still known. Approximate al- 
gorithms such as the i?-algorithm j|] have been 
used in these cases. 

Several exact algorithms are proposed for two- 
or single-flavor dynamical KS fermions In 
this paper, we further study the idea by Horvath 
et al. in the case of the polynomial Hybrid 
Monte Carlo (PHMC) algorithm. We develop two 
types of the PHMC algorithm depending on the 



choice of the effective action derived through the 
polynomial approximation. We compare the com- 
putational cost of these two PHMC algorithms. 
We investigate the property of the algorithm on 
several lattice sizes in Nf=2 case. We found that 
our algorithm shows satisfactory efficiency on a 
le'' lattice with /3=5.7, to=0.02 {mps/mv^Om). 

2. Algorithm 

We construct two types of the PHMC algo- 
rithm, which are refereed to as case A and case 
B. Introducing a polynomial approximation and 
pseudo-fermion field, the partition function can 
be generally rewritten in the following form: 



Z=J VUV(t)lV<j)o det[VF(^)[i)o 
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y^f,-sAu]-si''>lu,'Pl<Po]^ (2) 

where Sg is a lattice gauge action, 0o is pseudo- 
fermion field living only on odd sites. Sq and 
are the pseudo-fermion action and the cor- 
rection matrix respectively. The superscript {X) 
takes (A) or (B) depending on the type of the 
PHMC algorithm as follows. 
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Case A: We approximate l), 



-Nf/8 



by a N^^) 

•J poly 



order polynomial Pj^(a) [Doo] ■ The pseudo- 

poly 

fermion action and the correction term become 

poly 

W^^^[Doo] = Doo{P^iA,[Doo]f^'''. (2) 



Case B: We approximate l), 

r(B) 
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by an even- 



order N^l^J polynomial Pj^(b)[Doo]- The pscudo- 

poly 

fermion action and the correction term can be 
written as 



|g^(B,[^oo]0ol', 

DooiP^iB.iDoo])'/"'- , (3) 



where Qj^{b) is the order polynomial de- 

fined by Pj^(b)[Doo]^\Q j^(b)[Doo]\'^ ■ 

poly poly 

The KS-fermion operator Dqo is even-odd pre- 
conditioned as 000 = ^00— Moo with A^=2Ainax/ 

(4to2+2A2„^^) and Moo=2M<,eMeo/A^a^-hloo. 
Amax is chosen so that the all eigenvalues of Moo 
fall into the region [—1, 1]. Moe {Mf,o) is the usual 
KS hopping matrix from even (odd) to odd (even) 
sites. 

For both cases, the algorithm takes the follow- 
ing two steps; (i) perform the HMC algorithm 
according to the effective action Eq. (|^) or (H), 
(ii) when the HMC Metropolis test is accepted, 
apply the Kennedy-Kuti noisy Metropolis test to 
incorporate the correction term W^'^'>[Doo]]^'^'^ ■ 
Thus we obtain two types of the PHMC algorithm 
depending on the choice of the effective action. 

The acceptance probability of the noisy- 
Metropolis test is defined by 

Pnmp[(7^C/'] =min[l,e-'*^], 
dS^ClW^'''>[D'oo]-''^/'Co-\Vo\', (4) 

where (o=W'^^'^[Doo]^f^^^Vo with a Gaussian 
noise vector rjo. W^'^''[Doo] is calculated on an 
initial configuration and W^-^^D'oo] is on a trial 
configuration generated by the preceding HMC 
algorithm. 

The fractional power of the correction matrix 
W'^^^ is taken by the Lanczos based Krylov sub- 



space method proposed by Borigi [|6|. We modi- 
fied his algorithm suitable to our purpose. As in- 
dicated by Borigi we employ CG based stopping 
criterion for the Lanczos based method. 

3. Cost estimate 

The computational cost is counted as the num- 
ber of multiplication of the hopping matrix to 
evolve the algorithm unit trajectory. We employ 
single leapfrog integration scheme for the molec- 
ular dynamics (MD) step. We roughly estimate 
it as 



- (2<.l - 1) X <^ 
+3 X ((8/A^/) X A^(^) + 1) X Ar(^\ 

CostB = {n12 1) X iVlS 

+3x((4/A^^.)xA^S + l)xArg), 



CostA 



where A'^md is the number of MD step and A^qq 
the number of iteration of CG algorithm. The CG 
algorithm is used to generate (j)o with the global 
heat-bath method and in the Lanczos based algo- 
rithm for the noisy Metropolis test. 

Now we compare the costs by specifying N^^^^ 

and Np^i'y. For this purpose we employ the 
Chebyshev polynomial approximation; 

x-^ = (1 + X'y)-^ ^ Pn^^Jx] = ^fc^feM' (5) 



where Tk is the fc-th order Chebyshev polynomial, 
x=Doo and y=—Moo can be read of. The coeffi- 
cients Cfc are calculated as usual. 

Figure |l| shows the cost N^oiy/s depen- 
dence of the integrated residual defined by 

^=U-, dy\x{PN,^,M)"' - IH'/' for each s. 
We observe that as decreasing s by factor 1/2 the 
cost Npoiy/s increases by factor 2 at a constant 
This is nothing but Npoiy does not depend 
on the choice of s and we obtain nI^), -A^lf!! at 

poiy poiy 

a constant approximation level. 

Using this relation and assuming N^=N^ 

and Njj^2=Ni^B^ we find iVcostA'-2A^costB. We 
employ the Chebyshev polynomial and the pla- 
quette gauge action, and apply the case B PHMC 
algorithm for numerical simulations. 
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Figure 1. Npoiy/s dependence of the integrated 
residual ^/W with 1 - = 1/1000. 

Table 1 

Numerical results on a 16"* lattice at /3=5.7 and 
m=0.02. Amax= 2.28 is employed. 



JVpoly 

[dt, ATmd] 


300 


400 


500 


[0.02,50] 


[0.02,50] 


[0.02, 50] 


Traj. 


1700 


1050 


800 


<^'> 


0.577099(46) 


0.577130(46) 


0.577023(43) 



4. Results 

Figure || shows the MD step size dependence 
of the averaged plaquette on a 8'^ x 4 lattice 
at /3=5.26 and m=0.025. The results with the 
PHJVIC algorithm of iVpoj^ =200 do not depend 
on dt as it should be and produce the consistent 
result to that in the zero MD step size limit of 
the i?-algorithm. Although we do not show the 
results on the Npoiy dependence, the results are 
independent of the choice of Npoiy . 

In Table |l| we show the numerical results on a 
16'' lattice at (3=5.7 and to=0.02. The averaged 
plaquette are independent of Npdy and consistent 
with each other as expected. On the other hand, 
the i?-algorithm yields (F) =0.577261 (49) 0, 
which differs from ours by ~2ct. This indicates a 
potential systematic error for the i?-algorithm at 
finite MD step size. The computational time for 
A^po;y=300 with m=0.02 (which corresponds to 
mps/TOv~0.69 0) was measured as 112 sec. to 
achieve unit trajectory with 14 GFlops sustained 
speed of SR8000 at KEK. 

Consequently we conclude that the PHMC al- 
gorithm we constructed works on a moderately 
large lattice size with rather heavy quark masses 
with reasonable computational cost. The algo- 
rithm also works with a single-flavor fermion. We 
emphasize that because the PHMC algorithm is 
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Figure 2. MD step size dt=l/Nyu:) dependence 
of the averaged plaquette (P) on the small size 
lattice. Npoiy=2Q{), iVMD=25, and Amax=2.37 are 
employed for the PHMC. 

exact one, it must be a promising algorithm for 
future realistic simulations. 
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